Descriptive microbiome stats and cleaner code for dissertation / publication.
library(magrittr)
library(phyloseq)
library(fs)
library(tidyverse)
library(vegan)
library(rstatix)
library(ggpubr)
library(microViz)
library(pals)
library(RColorBrewer)
library(colorBlindness)
library(microbiome)
library(ggrepel)
library(DT)
library(plotly)
library(ape)
library(picante)
library(phytools)
library(DESeq2)
library(parallel)
library(doBy)
library(ragg)
library(cowplot)
rm(list=ls()) # clear environment
git.dir = file.path("/hpc/group/dmc/Eva/GLP_KO_Microbiome_final") # CHANGE ME
ps.rds <- file.path(git.dir, "ps.GLP1.rds")
ps <- read_rds(ps.rds)
ps
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 234 samples ]
## sample_data() Sample Data: [ 234 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]
fig.dir = file.path(git.dir, "figures")
map.file = file.path(git.dir, "Metadata_GLP1.csv")
meta.df = read_csv(map.file, show_col_types = FALSE)
Note that this phyloseq object is already filtered such that only samples with reads > 1000 are available. See here:
sample_sums(ps) %>% as.data.frame() %>%
dplyr::rename(Reads = ".") %>%
arrange(Reads)
The names are hard to interpret because they are the SRA accession numbers. Let’s change the names in the phyloseq object.
phyloseq::sample_names(ps) <- phyloseq::sample_data(ps)$Label
meta.df$Genotype <- factor(meta.df$Genotype, levels = c("WT", "KO"))
meta.df$Sex <- factor(meta.df$Sex, levels = c("Female", "Male"))
meta.df$Diet <- factor(meta.df$Diet, levels = c("Normal_Chow", "High_Fat_Diet"))
meta.df$Intranasal_Treatment <- factor(meta.df$Intranasal_Treatment, levels = c("PBS", "HDM"))
meta.df$Surgery <- factor(meta.df$Surgery, levels = c("None", "Sham", "VSG"))
meta.df$Group <- factor(meta.df$Group, levels = c("NA", "Control", "1", "2"))
meta.df %>% dplyr::filter(Label %in% rownames(phyloseq::sample_data(ps))) %>%
mutate(Label_graph = Label) %>%
column_to_rownames("Label") -> meta.ps
ps
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 234 samples ]
## sample_data() Sample Data: [ 234 samples by 42 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]
sample_data(ps) <- meta.ps
ps
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 234 samples ]
## sample_data() Sample Data: [ 234 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]
Genotype, Sex, Diet,
Intranasal_Treatment, and Surgery are
appropriately factors now.
# Only include samples
subset_samples(ps, Type == "True Sample") -> ps_sam
# remove samples with potential labeling issues
ps_sam %>% subset_samples(row.names(ps_sam@sam_data) != "294" & row.names(ps_sam@sam_data) != "30") -> ps_sam
# Subset by sample type
ps_sam %>% subset_samples(Sample_type == "Lung Tissue") -> ps.LT.asv
ps_sam %>% subset_samples(Sample_type == "Fecal") -> ps.F.asv
ps_sam
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 221 samples ]
## sample_data() Sample Data: [ 221 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]
ps.LT.asv
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 44 samples ]
## sample_data() Sample Data: [ 44 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]
ps.F.asv
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 177 samples ]
## sample_data() Sample Data: [ 177 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]
# Relative abundance
ps_sam.ts <- transform_sample_counts(ps_sam, function(x) x/sum(x))
ps.LT.asv.ts <- transform_sample_counts(ps.LT.asv, function(x) x/sum(x))
ps.F.asv.ts <- transform_sample_counts(ps.F.asv, function(x) x/sum(x))
Let’s subset the metadata to those that are in the
ps_sam object.
meta.df %>%
filter(meta.df$Label %in% rownames(ps_sam@sam_data)) -> meta.sam
meta.sam %>% filter(is.na(Genotype) == TRUE & is.na(Intranasal_Treatment) == TRUE & is.na(Surgery) == TRUE)
rarecurve(as.data.frame(otu_table(ps_sam)), step = 100, cex = 0.5, tidy = TRUE) -> rare_ps #put in ggplot2 friendly dataframe
## Warning in rarecurve(as.data.frame(otu_table(ps_sam)), step = 100, cex = 0.5, :
## most observed count data have counts 1, but smallest count is 2
meta.sam %>% dplyr::select(Label, Sample_type, Mouse, Week, DNA_ng_ul_raw) %>%
dplyr::rename(Site = Label) %>%
right_join(rare_ps, by = "Site") -> rare_ps
ggplot(data=rare_ps, aes(x=Sample, y=Species, group=Site)) +
geom_line() + theme_bw() + labs(x = "Number of Reads", y= "Species Richness") +
facet_wrap(~Sample_type, scales = "free") -> p
p
ggplotly(p, tooltip = "Site")
ps_sam %>% tax_glom("Phylum") %>%
transform_sample_counts(function(x) x/sum(x)) -> ps.phylum.ts
taxa_names(ps.phylum.ts) <- tax_table(ps.phylum.ts)[, 'Phylum']
ps.phylum.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 20 taxa and 221 samples ]
## sample_data() Sample Data: [ 221 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 20 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 20 tips and 19 internal nodes ]
## refseq() DNAStringSet: [ 20 reference sequences ]
plot_bar(ps.phylum.ts,x = "Label_graph", fill="Phylum") +
geom_bar(aes(color=Phylum, fill=Phylum), stat="identity", position="stack") +
theme_bw() +
labs(y="Relative Abundance", x = "Sample") +
coord_cartesian(ylim = c(0,1), expand=0) +
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) +
theme(axis.text.x=element_blank()) +
facet_wrap(~Sample_type, scales = "free") -> p
p
ps_melt(ps.phylum.ts) -> phylum.ts.melt
phylum.ts.melt %>%
filter(Sample_type == "Lung Tissue") %>%
doBy::summary_by(Abundance ~ OTU, FUN = median) %>%
arrange(desc(Abundance.median))
phylum.ts.melt %>%
filter(Sample_type == "Fecal") %>%
doBy::summary_by(Abundance ~ OTU, FUN = median) %>%
arrange(desc(Abundance.median))
ps.phylum.ts %>%
subset_samples(Sample_type == "Fecal") %>%
ps_arrange(desc(Firmicutes), desc(Bacteroidota), .target = "otu_table") -> ps.fecal.phylum.ts
sample_data(ps.fecal.phylum.ts)$Sample_Order <- c(1:nrow(sample_data(ps.fecal.phylum.ts)))
ps.fecal.phylum.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 20 taxa and 177 samples ]
## sample_data() Sample Data: [ 177 samples by 41 sample variables ]
## tax_table() Taxonomy Table: [ 20 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 20 tips and 19 internal nodes ]
## refseq() DNAStringSet: [ 20 reference sequences ]
ps.phylum.ts %>%
subset_samples(Sample_type == "Lung Tissue") %>%
ps_arrange(desc(Proteobacteria), desc(Actinobacteriota), .target = "otu_table") -> ps.lt.phylum.ts
sample_data(ps.lt.phylum.ts)$Sample_Order <- c(1:nrow(sample_data(ps.lt.phylum.ts)))
ps.lt.phylum.ts
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 20 taxa and 44 samples ]
## sample_data() Sample Data: [ 44 samples by 41 sample variables ]
## tax_table() Taxonomy Table: [ 20 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 20 tips and 19 internal nodes ]
## refseq() DNAStringSet: [ 20 reference sequences ]
ps_melt(ps.fecal.phylum.ts) %>%
ggplot(aes(fill=Phylum, y=Abundance, x=Sample_Order)) +
geom_bar(position="fill", stat="identity", width = 1) +
theme_bw() +
labs(y="Relative Abundance", x = "Sample") +
coord_cartesian(ylim = c(0,1), expand=0) +
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) +
facet_wrap(~Sample_type) +
theme(axis.text.x=element_blank()) -> p.fecal
p.fecal
ps_melt(ps.lt.phylum.ts) %>%
ggplot(aes(fill=Phylum, y=Abundance, x=Sample_Order)) +
geom_bar(position="fill", stat="identity", width = 1) +
theme_bw() +
labs(y="Relative Abundance", x = "Sample") +
coord_cartesian(ylim = c(0,1), expand=0) +
scale_color_manual(values=SteppedSequential5Steps) +
scale_fill_manual(values=SteppedSequential5Steps) +
facet_wrap(~Sample_type) +
theme(axis.text.x=element_blank()) -> p.lt
p.lt
ragg::agg_jpeg(file.path(fig.dir, "fig_RelAb_fecal.jpeg"), width = 10, height = 6, units = "in", res = 600)
p.fecal
dev.off()
## png
## 2
ragg::agg_jpeg(file.path(fig.dir, "fig_RelAb_lt.jpeg"), width = 5, height = 6, units = "in", res = 600)
p.lt
dev.off()
## png
## 2
ps.LT.asv %>% estimate_richness(measures = c("Observed", "Shannon")) -> LT_richness
LT_richness$Faith <- pd(otu_table(ps.LT.asv), phy_tree(ps.LT.asv), include.root = FALSE)$PD
LT_richness %>% rownames_to_column(var = "Label") -> LT_richness
gsub("\\.", "-", LT_richness$Label) -> LT_richness$Label
meta.df %>%
filter(meta.df$Label %in% LT_richness$Label) %>% # merge with metadata
select(c(Label:DNA_ng_ul_raw, "Group")) %>%
right_join(LT_richness, by = "Label") -> LT_richness
LT_richness %>% pivot_longer(
cols = contains(c("Observed", "Shannon", "Faith")),
names_to = c("Measurement")
) -> LT_richness_stat # dataframe that works well for ggplot2() and anova_test()
LT_richness_stat$Mouse <- factor(LT_richness_stat$Mouse)
LT_richness_stat$Measurement <- factor(LT_richness_stat$Measurement)
LT_richness_stat$Cohort <- factor(LT_richness_stat$Cohort)
LT_richness %>% nrow()
## [1] 44
LT_richness_stat %>% nrow()
## [1] 132
44*3
## [1] 132
We are most interested in intranasal treatment, diet, and genotype if possible. We also note that there were some mice whose lung tissues were collected on week 10 vs. week 13.
LT_richness %>% filter(Surgery == "None") -> lt.richness.nosurgery
lt.richness.nosurgery %>%
dplyr::count(Week, Genotype, Intranasal_Treatment) %>% arrange(n)
lt.richness.nosurgery %>%
dplyr::count(Week, Intranasal_Treatment, Diet) %>% arrange(n)
Let’s graph:
LT_richness_stat %>% filter(Surgery == "None") -> lt.richness.nosurgery.long
lt.richness.nosurgery.long %>%
ggplot(aes(x = Intranasal_Treatment, y = value)) +
geom_boxplot() + geom_point() +
facet_wrap(vars(Measurement), scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = Genotype, y = value)) +
geom_boxplot() + geom_point() +
facet_wrap(vars(Measurement), scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = Diet, y = value)) +
geom_boxplot() + geom_point() +
facet_wrap(vars(Measurement), scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = as.character(Week), y = value)) +
geom_boxplot() + geom_point() +
facet_wrap(vars(Measurement), scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = Genotype, y = value)) +
geom_boxplot() + geom_point() +
facet_grid(vars(Measurement), vars(Intranasal_Treatment),
scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = Diet, y = value)) +
geom_boxplot() + geom_point() +
facet_grid(vars(Measurement), vars(Intranasal_Treatment),
scales = "free") + theme_bw()
lt.richness.nosurgery.long %>%
ggplot(aes(x = as.character(Week), y = value)) +
geom_boxplot() + geom_point() +
facet_grid(vars(Measurement), vars(Intranasal_Treatment),
scales = "free") + theme_bw()
Graphically, it looks like intranasal treatment, week, and genotype may be of interest.
lt.richness.nosurgery %>% nrow()
## [1] 26
lt.richness.nosurgery %>%
dplyr::count(Week, Intranasal_Treatment) %>% arrange(n)
lt.richness.nosurgery.long %>%
group_by(Week, Intranasal_Treatment, Measurement) %>%
shapiro_test(value) %>%
add_significance()
Mainly indistinguishable from normal distribution.
lt.richness.nosurgery.long %>%
group_by(Week, Intranasal_Treatment, Measurement) %>%
identify_outliers(value) -> lt.nosurgery.outliers
lt.nosurgery.outliers
Outliers only in week 10 groups.
meta.sam %>% filter(Sample_type == "Lung Tissue") %>%
filter(Label %in% lt.nosurgery.outliers$Label) %>%
dplyr::select(Label, Sample_type:DNA_ng_ul_raw, Notes)
Interesting, maybe genotype or sex is playing a role here?
lt.richness.nosurgery.long %>%
ggqqplot(x = "value") + facet_grid(vars(Measurement), vars(Intranasal_Treatment), scales = "free")
lt.richness.nosurgery.long %>%
group_by(Measurement) %>%
anova_test(
dv = value,
between = c(Week, Intranasal_Treatment, Genotype),
) %>% get_anova_table() %>%
adjust_pvalue(method = "holm") %>%
add_significance()
Nothing is significant, sadly. Let’s graph with intranasal treatment, as that had the greatest effect on host respiratory health metrics.
lt.richness.nosurgery.long %>%
group_by(Measurement) %>%
t_test(value ~ Intranasal_Treatment) %>%
adjust_pvalue(method = "holm") %>%
add_significance() -> lt.richness.nosurgery.res
lt.richness.nosurgery.res
# Change facet label names
alpha.labs <- c("Faith's PD", "Observed Richness", "Shannon")
names(alpha.labs) <- c("Faith", "Observed", "Shannon")
#lt.richness.nosurgery.res$y.position <- c(40, 225, 4.7)
lt.richness.nosurgery.long %>%
ggplot(aes(x = Intranasal_Treatment, y = value, color = Intranasal_Treatment))+
geom_boxplot() +
geom_point() +
facet_wrap(~Measurement, scales = "free_y",
labeller = labeller(Measurement = alpha.labs)) + expand_limits(y = 0) +
theme_bw(base_size = 14) + scale_fill_brewer(palette="RdBu") +
labs(y="Alpha Diversity Value", x="Intranasal Treatment", color = "Intranasal Treatment") +
theme(legend.position = "top") -> lt.alpha.group1.week
# stat_pvalue_manual(lt.richness.nosurgery.res, label = "p.adj.signif") -> lt.alpha.group1.week
lt.alpha.group1.week
For the after surgery timepoint for lung tissue, look at week 13 only and HFD.
LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13) %>%
dplyr::count(Surgery, Intranasal_Treatment)
LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13) %>%
dplyr::count(Surgery)
LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13) %>%
dplyr::count(Intranasal_Treatment)
LT_richness %>% filter(Diet == "High_Fat_Diet" & Week == 13)%>% nrow()
## [1] 23
LT_richness_stat %>% filter(Diet == "High_Fat_Diet" & Week == 13) -> lt.richness.hfd
lt.richness.hfd %>%
group_by(Measurement, Surgery) %>%
shapiro_test(value) %>%
add_significance()
lt.richness.hfd %>% ggqqplot(x = "value") + facet_grid(vars(Measurement), vars(Surgery), scales = "free")
It looks like this is due to an outlier, but rarefaction curves looked like they plateaued. Proceed with caution.
lt.richness.hfd %>%
group_by(Measurement, Surgery) %>%
shapiro_test(value) %>%
add_significance()
lt.richness.hfd %>%
group_by(Measurement, Surgery) %>%
identify_outliers(value)
lt.richness.hfd %>%
ggplot(aes(x = Surgery, y = value)) +
geom_boxplot() + geom_point() +
facet_wrap(vars(Measurement), scales = "free") + theme_bw()
lt.richness.hfd %>%
group_by(Measurement) %>%
anova_test(
dv = value,
between = c(Surgery, Intranasal_Treatment),
) %>% get_anova_table() %>%
adjust_pvalue(method = "holm") %>%
add_significance()
Nothing is significant. Let’s make a graph showing surgery:
lt.richness.hfd %>%
ggplot(aes(x = Surgery, y = value, color = Surgery)) +
geom_boxplot() + geom_point() +
facet_wrap(~Measurement, scales = "free_y",
labeller = labeller(Measurement = alpha.labs)) + expand_limits(y = 0) +
theme_bw(base_size = 14) + scale_color_brewer(palette="Dark2") +
labs(y="Alpha Diversity Value", x="Surgery", color = "Surgery") +
theme(legend.position = "top") -> fig.lt.hfd
fig.lt.hfd
plot_grid(
lt.alpha.group1.week, fig.lt.hfd,
labels = "AUTO", ncol = 2, label_size = 20, scale = 0.95,
rel_widths = c(1, 1.3)
) -> fig.lt.alpha
fig.lt.alpha
ragg::agg_jpeg(file.path(fig.dir, "fig_lt_alpha.jpeg"), width = 14, height = 5, units = "in", res = 600, scaling = 1)
fig.lt.alpha
dev.off()
## png
## 2
For PERMANOVA, let’s have n>=5 for each group.
ps.LT.asv.ts %>% subset_samples(Surgery == "None") -> ps.lt.ts.nosurgery
# Calculate ditances
ps.lt.ts.nosurgery %>% phyloseq::distance(method = "unifrac") -> LT_nosurgery_uni
ps.lt.ts.nosurgery %>% phyloseq::distance(method = "wunifrac") -> LT_nosurgery_weighted_uni
# Ordinate distances
LT_nosurgery.uni_ord <- ordinate(ps.lt.ts.nosurgery, method = "PCoA", distance = "unifrac")
LT_nosurgery.weighted_uni_ord <- ordinate(ps.lt.ts.nosurgery, method = "PCoA", distance = "wunifrac")
# Extract dataframe
ps.lt.ts.nosurgery@sam_data %>% data.frame() -> LT_nosurgery.df
LT_nosurgery.df %>% dplyr::count(Genotype, Intranasal_Treatment) %>% arrange(n)
# Set seed for reproducibility
# UniFrac
set.seed(123)
adonis2(formula = LT_nosurgery_uni ~ Intranasal_Treatment * Genotype, data = LT_nosurgery.df, permutations = 9999, parallel = getOption("mc.cores"))
# Weighted UniFrac
set.seed(123)
adonis2(formula = LT_nosurgery_weighted_uni ~ Intranasal_Treatment * Genotype, data = LT_nosurgery.df, permutations = 9999, parallel = getOption("mc.cores"))
Not significant.
plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.uni_ord, type = "samples", color = "Intranasal_Treatment") + theme_bw() +
labs(title = "UniFrac") +
scale_color_brewer(palette="Dark2")
plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.weighted_uni_ord, type = "samples", color = "Intranasal_Treatment") + theme_bw() +
labs(title = "Weighted UniFrac") +
scale_color_brewer(palette="Dark2")
plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.uni_ord, type = "samples", color = "Intranasal_Treatment") + theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") +
theme(legend.position = "top") -> p.lt.nosurgery.uni
p.lt.nosurgery.uni
plot_ordination(ps.lt.ts.nosurgery, LT_nosurgery.weighted_uni_ord, type = "samples", color = "Intranasal_Treatment") +
theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") +
theme(legend.position = "top") -> p.lt.nosurgery.wuni
p.lt.nosurgery.wuni
plot_grid(p.lt.nosurgery.uni, p.lt.nosurgery.wuni,
ncol = 2, labels = "AUTO", label_size = 20) -> p.lt.nosurgery
p.lt.nosurgery
ragg::agg_jpeg(file.path(fig.dir, "fig_lt_nosurgery_beta.jpeg"), width = 10, height = 5, units = "in", res = 600, scaling = 1.1)
p.lt.nosurgery
dev.off()
## png
## 2
Only look at at HFD at week 13.
ps.LT.asv.ts %>% subset_samples(Diet == "High_Fat_Diet") %>% subset_samples(Week == "13") -> ps.lt.ts.hfd
# Calculate ditances
ps.lt.ts.hfd %>% phyloseq::distance(method = "unifrac") -> LT_hfd_uni
ps.lt.ts.hfd %>% phyloseq::distance(method = "wunifrac") -> LT_hfd_weighted_uni
# Ordinate distances
LT_hfd.uni_ord <- ordinate(ps.lt.ts.hfd, method = "PCoA", distance = "unifrac")
LT_hfd.weighted_uni_ord <- ordinate(ps.lt.ts.hfd, method = "PCoA", distance = "wunifrac")
# Extract dataframe
ps.lt.ts.hfd@sam_data %>% data.frame() -> LT_hfd.df
LT_hfd.df %>% dplyr::count(Surgery) %>% arrange(n)
# Set seed for reproducibility
# UniFrac
set.seed(123)
adonis2(formula = LT_hfd_uni ~ Surgery, data = LT_hfd.df, permutations = 9999, parallel = getOption("mc.cores"))
# Weighted UniFrac
set.seed(123)
adonis2(formula = LT_hfd_weighted_uni ~ Surgery, data = LT_hfd.df, permutations = 9999, parallel = getOption("mc.cores"))
NS.
plot_ordination(ps.lt.ts.hfd, LT_hfd.uni_ord, type = "samples", color = "Surgery") + theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") +
theme(legend.position = "top") -> p.lt.hfd.uni
p.lt.hfd.uni
plot_ordination(ps.lt.ts.hfd, LT_hfd.weighted_uni_ord, type = "samples", color = "Surgery") +
theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") +
theme(legend.position = "top") -> p.lt.hfd.wuni
p.lt.hfd.wuni
plot_grid(p.lt.hfd.uni, p.lt.hfd.wuni,
ncol = 2, labels = "AUTO", label_size = 20) -> p.lt.hfd
p.lt.hfd
ragg::agg_jpeg(file.path(fig.dir, "fig_lt_hfd_beta.jpeg"), width = 10, height = 5, units = "in", res = 600, scaling = 1.1)
p.lt.hfd
dev.off()
## png
## 2
ps.F.asv
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 177 samples ]
## sample_data() Sample Data: [ 177 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]
ps.F.asv %>% estimate_richness(measures = c("Observed", "Shannon")) -> F_richness
## Warning in estimate_richness(., measures = c("Observed", "Shannon")): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.
F_richness$Faith <- pd(otu_table(ps.F.asv), phy_tree(ps.F.asv), include.root = FALSE)$PD
F_richness %>% rownames_to_column(var = "Label") -> F_richness
gsub("X", "", F_richness$Label) -> F_richness$Label # remove X in front of sample names
meta.df %>%
filter(meta.df$Label %in% F_richness$Label) %>% # merge with metadata
select(c(Label, Genotype:Surgery, Kit:Feces_weight, Week, Cohort, Mouse, Notes)) %>%
right_join(F_richness, by = "Label") -> F_richness
F_richness %>% group_by(Mouse, Week) %>% filter(n()>1) -> F_richness.dup
F_richness.dup
# calculate averages for duplicates
F_richness.dup %>%
dplyr::select(Label, Mouse, Week, Observed:Faith) %>% group_by(Mouse) %>%
dplyr::summarise(Observed = mean(Observed),
Shannon = mean(Shannon),
Faith = mean(Faith)) -> F_richness.dup.calc
# add week information for duplicates
F_richness.dup %>% dplyr::select(Mouse, Week) %>% unique() %>%
ungroup() %>%
right_join(F_richness.dup.calc, by = "Mouse") -> F_richness.dup.calc
# add metadata for duplicates
meta.df %>% filter(Sample_type == "Fecal") %>%
filter(Mouse %in% F_richness.dup.calc$Mouse) %>%
select(c(Genotype:Surgery, Week, Cohort, Mouse, Notes)) %>%
right_join(F_richness.dup.calc, by = c("Mouse", "Week")) %>%
unique() -> F_richness.dup.calc
F_richness.dup.calc
# Clean up: remove duplicates from original dataframe, then add calculated values.
F_richness %>% filter(!Label %in% F_richness.dup$Label) %>%
bind_rows(F_richness.dup.calc) -> F_richness
F_richness %>% filter(is.na(Diet) == FALSE) -> F_richness
F_richness %>% pivot_longer(
cols = contains(c("Observed", "Shannon", "Faith")),
names_to = c("Measurement")
) -> F_richness_stat # dataframe that works well for ggplot2() and anova_test()
F_richness_stat$Mouse <- factor(F_richness_stat$Mouse)
F_richness_stat$Measurement <- factor(F_richness_stat$Measurement)
F_richness_stat$Cohort <- factor(F_richness_stat$Cohort)
I want to see if high fat diet resulted in a significant decline in
alpha diversity, as has been reported before. To do this, I will compare
wk 0 and wk 10 alpha diversity for all mice
based on diet (to increase n). I also want to see if genotype
and sex had an effect.
F_richness %>% filter(Week != 13) -> F_richness.nosurgery
F_richness.nosurgery%>%
dplyr::count(Week, Diet, Genotype, Sex) %>% arrange(n)
F_richness.nosurgery %>% distinct(Mouse) %>% nrow()
## [1] 69
F_richness.nosurgery %>% nrow()
## [1] 103
F_richness_stat %>% filter(Week != 13) -> F_richness.nosurgery.long
F_richness.nosurgery.long %>%
group_by(Measurement, Week, Diet, Genotype, Sex) %>%
shapiro_test(value) %>%
add_significance()
Mainly indistinguishable from normal distribution. Let’s graph the two that are distinguishable from normal distribution:
F_richness.nosurgery.long %>%
filter(Sex == "Female" & Genotype == "WT" & Diet == "Normal_Chow") %>%
ggqqplot(x = "value") + theme_bw() + facet_grid(vars(Measurement), vars(Week), scales = "free")
It looks like it’s driven by outliers, but otherwise looks quite good.
F_richness.nosurgery.long %>%
group_by(Measurement, Week, Diet, Genotype, Sex) %>%
identify_outliers(value)
The feces weight, gDNA concentrations, and notes (or the lack thereof) for even extreme outliers look good.
F_richness.nosurgery.long %>%
group_by(Measurement) %>%
anova_test(dv = value,
wid = Mouse,
within = Week,
between = c(Diet, Genotype, Sex)) %>%
get_anova_table() %>%
as.data.frame() %>%
adjust_pvalue(method = "holm") %>%
add_significance() -> F.0.10.res
F.0.10.res
F.0.10.res %>% filter(p.adj < 0.05)
Only diet is significant. Let’s redo analysis and remove genotype and sex as factors, as they were not significant.
F_richness.nosurgery%>%
dplyr::count(Week, Diet) %>% arrange(n)
F_richness.nosurgery %>% distinct(Mouse) %>% nrow()
## [1] 69
F_richness.nosurgery %>% nrow()
## [1] 103
F_richness.nosurgery.long %>%
group_by(Measurement, Week, Diet) %>%
shapiro_test(value) %>%
add_significance()
F_richness.nosurgery.long %>%
group_by(Measurement, Week, Diet) %>%
identify_outliers(value)
None are extreme outliers. Proceed:
F_richness.nosurgery.long %>%
group_by(Measurement) %>%
anova_test(dv = value,
wid = Mouse,
within = Week,
between = c(Diet)) %>%
get_anova_table() %>%
as.data.frame() %>%
adjust_pvalue(method = "holm") %>%
add_significance() -> F.0.10.res.diet
F.0.10.res.diet
F.0.10.res.diet %>% filter(p.adj < 0.05)
Let’s break down the diet:week interaction in Faith and Shannon diversity metrics. First, let’s graph what they look like:
F_richness.nosurgery.long$Week <- factor(F_richness.nosurgery.long$Week, levels = c("0", "10"))
# Change facet label names
alpha.labs <- c("Faith's PD", "Observed Richness", "Shannon")
names(alpha.labs) <- c("Faith", "Observed", "Shannon")
ggplot(F_richness.nosurgery.long, aes(x = Week, y = value, color = Diet))+
geom_boxplot(position=position_dodge(width = 0.8)) + geom_point(position=position_dodge(width = 0.8)) +
facet_wrap(~Measurement, scales = "free", labeller = labeller(Measurement = alpha.labs)) +
labs(x="Week", y="Alpha Diversity Metric Value", fill="Diet") + expand_limits(y = 0) + theme_bw() +
theme(legend.position = "top")
# other way to test / graph
ggplot(F_richness.nosurgery.long, aes(x = Diet, y = value, color = Week))+
geom_boxplot(position=position_dodge(width = 0.8)) + geom_point(position=position_dodge(width = 0.8)) +
facet_wrap(~Measurement, scales = "free", labeller = labeller(Measurement = alpha.labs)) +
labs(x="Diet", y="Alpha Diversity Metric Value", fill="Week") + expand_limits(y = 0) + theme_bw() +
theme(legend.position = "top")
Since Diet:Week wasn’t significant for Observed, let’s correct using Holm adjustment.
F_richness.nosurgery.long %>% # raw values only
group_by(Measurement, Week) %>%
t_test(value ~ Diet) %>%
adjust_pvalue(method = "holm") %>%
add_significance() -> F_wk0_10_ttest
F_wk0_10_ttest
F_wk0_10_ttest %>%
dplyr::select(Measurement, Week, n1:n2, p.adj, p.adj.signif)
# plot mean +/- standard error
fun_se <- function(x){
c(std_error = sd(x)/sqrt(length(x)), mean = mean(x))
}
# calculate mean and se for each week, measurement, and diet
F_richness.nosurgery.long %>%
doBy::summary_by(value ~ Measurement + Week + Diet, FUN = fun_se) -> F_richness.wk0.10.graph
ggplot(F_richness.wk0.10.graph, aes(x = Week, y = value.mean, color = Diet))+
geom_line(aes(group=Diet)) +
geom_errorbar(aes(ymin=value.mean-value.std_error, ymax=value.mean+value.std_error, width=0.2)) +
theme_bw(base_size = 14) + scale_fill_brewer(palette="RdBu") +
facet_wrap(~Measurement, scales = "free", labeller = labeller(Measurement = alpha.labs)) +
labs(x="Week", y="Alpha Diversity Metric Value", fill="Diet") + expand_limits(y = 0) +
theme(legend.position = "top") -> fig.alpha.0.10
fig.alpha.0.10
ragg::agg_jpeg(file.path(fig.dir, "fig_alpha_0to10.jpeg"), width = 7, height = 5, units = "in", res = 600, scaling = 1.2)
fig.alpha.0.10
dev.off()
## png
## 2
F_richness %>%
filter(Diet == "High_Fat_Diet" & Week == 13) -> f.richness.hfd
f.richness.hfd %>% dplyr::count(Surgery) %>% arrange(n)
F_richness %>%
filter(Diet == "High_Fat_Diet" & Week == 13) %>% distinct(Mouse) %>% nrow()
## [1] 45
f.richness.hfd %>% dplyr::count(Surgery, Genotype) %>% arrange(n)
f.richness.hfd %>% dplyr::count(Surgery, Genotype, Sex) %>% arrange(n)
F_richness_stat %>%
filter(Diet == "High_Fat_Diet" & Week == 13) -> f.richness.hfd
f.richness.hfd %>%
group_by(Measurement, Surgery, Genotype) %>%
shapiro_test(value) %>%
add_significance()
f.richness.hfd %>%
group_by(Measurement, Surgery, Genotype) %>%
identify_outliers(value)
f.richness.hfd %>% filter(Genotype == "KO") %>%
ggqqplot(x = "value") + facet_grid(vars(Measurement), vars(Surgery), scales = "free")
The QQ plots look quite good, except for the few outliers.
f.richness.hfd %>%
group_by(Measurement) %>%
anova_test(dv = value,
wid = Mouse,
between = c(Surgery, Genotype, Sex)) %>%
get_anova_table() %>%
as.data.frame() %>%
adjust_pvalue(method = "holm") %>%
add_significance() -> F.hfd.res
F.hfd.res
All ns. Let’s graph:
f.richness.hfd %>%
ggplot(aes(x=Surgery, y=value, color=Surgery)) +
geom_boxplot() +
geom_point() +
scale_fill_brewer(palette="RdBu") +
labs( y="Alpha Diversity Metric value", color="Surgery", x = "Surgery") +
facet_wrap(vars(Measurement), scales = "free", labeller = labeller(Measurement = alpha.labs)) +
theme_bw(base_size = 14) + expand_limits(y = 0) + theme(legend.position = "top") -> fig.alpha.hfd
fig.alpha.hfd
ragg::agg_jpeg(file.path(fig.dir, "fig_f_alpha_hfd.jpeg"), width = 7, height = 5, units = "in", res = 600)
fig.alpha.hfd
dev.off()
## png
## 2
ps.F.asv.ts %>% subset_samples(Surgery == "None") -> F_Diet
F_Diet %>%
phyloseq::distance(method = "unifrac") -> F_uni
F_Diet %>%
phyloseq::distance(method = "wunifrac") -> F_weighted_uni
F_Diet@sam_data %>% data.frame() -> F.df
F_uni_ord <- ordinate(F_Diet, method = "PCoA", distance = "unifrac")
F_weighted_uni_ord <- ordinate(F_Diet, method = "PCoA", distance = "wunifrac")
F.df %>% dplyr::count(Diet, Week) %>% arrange(n)
F.df %>% dplyr::count(Diet, Genotype, Week) %>% arrange(n)
F.df %>% nrow()
## [1] 120
F.df %>% distinct(Mouse) %>% nrow()
## [1] 58
# Set seed for reproducibility
set.seed(123)
adonis2(formula = F_uni ~ Diet * Genotype * Week, data = F.df, permutations = 9999, parallel = getOption("mc.cores")) -> F_diet_uni.res
F_diet_uni.res
set.seed(123)
adonis2(formula = F_weighted_uni ~ Diet * Genotype * Week, data = F.df, permutations = 9999, parallel = getOption("mc.cores")) -> F_diet_weighted_uni.res
F_diet_weighted_uni.res
As expected, diet is highly significant. Genotype is significant but not any of its interaction terms.
plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Diet") + theme_bw(base_size = 14) +
scale_fill_brewer(palette="RdBu") + theme(legend.position = "top") +
facet_wrap(~Week) + stat_ellipse() -> fig.diet.uni
plot_ordination(F_Diet, F_weighted_uni_ord, type = "samples", color = "Diet") +
theme_bw(base_size = 14) +
scale_fill_brewer(palette="RdBu") + theme(legend.position = "top", panel.spacing.x = unit(8, "mm")) +
facet_wrap(~Week) + stat_ellipse() -> fig.diet.w.uni
fig.diet.uni
fig.diet.w.uni
plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Genotype") + theme_bw(base_size = 14) +
scale_fill_brewer(palette="RdBu") + theme(legend.position = "top") + stat_ellipse() -> fig.geno.uni
plot_ordination(F_Diet, F_weighted_uni_ord, type = "samples", color = "Genotype") +
theme_bw(base_size = 14) +
scale_fill_brewer(palette="RdBu") + theme(legend.position = "top", panel.spacing.x = unit(8, "mm")) + stat_ellipse() -> fig.geno.w.uni
# UniFrac: diet*genotype was significant
plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Genotype") + theme_bw(base_size = 14) +
scale_fill_brewer(palette="RdBu") + facet_wrap(~Diet) +
theme(legend.position = "top") + stat_ellipse() -> fig.diet.geno.uni
fig.geno.uni
fig.diet.geno.uni
fig.geno.w.uni
plot_ordination(F_Diet, F_uni_ord, type = "samples", color = "Genotype") + theme_bw(base_size = 14) +
scale_fill_brewer(palette="RdBu") +
theme(legend.position = "top", panel.spacing.x = unit(5, "mm")) +
facet_grid(vars(Diet),vars(Week)) + stat_ellipse() -> fig.diet.geno.week.uni
fig.diet.geno.week.uni
plot_grid(
fig.alpha.0.10, fig.diet.uni,
fig.diet.w.uni, fig.diet.geno.week.uni,
labels = "AUTO", ncol = 2, label_size = 20
) -> fig.diet.summary
fig.diet.summary
ragg::agg_jpeg(file.path(fig.dir, "fig_diet_summary.jpeg"), width = 18, height = 12, units = "in", res = 600)
fig.diet.summary
dev.off()
## png
## 2
HFD: Effect of surgery; focus on wk 13 only.
ps.F.asv.ts %>% subset_samples(Diet == "High_Fat_Diet" & Week == 13) -> F_HFD
F_HFD
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 3304 taxa and 45 samples ]
## sample_data() Sample Data: [ 45 samples by 40 sample variables ]
## tax_table() Taxonomy Table: [ 3304 taxa by 6 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 3304 tips and 3303 internal nodes ]
## refseq() DNAStringSet: [ 3304 reference sequences ]
F_HFD %>%
phyloseq::distance(method = "unifrac") -> F_uni.HFD
F_HFD %>%
phyloseq::distance(method = "wunifrac") -> F_weighted_uni.HFD
F_HFD@sam_data %>% data.frame() -> F_HFD.df
F_uni_ord.HFD <- ordinate(F_HFD, method = "PCoA", distance = "unifrac")
F_weighted_uni_ord.HFD <- ordinate(F_HFD, method = "PCoA", distance = "wunifrac")
F_HFD.df %>% dplyr::count(Surgery) %>% arrange(n)
F_HFD.df %>% dplyr::count(Genotype) %>% arrange(n)
F_HFD.df %>% dplyr::count(Surgery, Genotype) %>% arrange(n)
F_HFD.df %>% nrow()
## [1] 45
# Set seed for reproducibility
set.seed(123)
adonis2(formula = F_uni.HFD ~ Surgery * Genotype, data = F_HFD.df, permutations = 9999, parallel = getOption("mc.cores"))
set.seed(123)
adonis2(formula = F_weighted_uni.HFD ~ Surgery * Genotype, data = F_HFD.df, permutations = 9999, parallel = getOption("mc.cores"))
Genotype is significant, surgery is not.
# stat_ellipse for the ones that are *
plot_ordination(F_HFD, F_uni_ord.HFD, type = "samples", color = "Genotype") + theme_bw(base_size = 14) +
theme(legend.position = "top") +
stat_ellipse() -> fig.hfd.uni.geno
fig.hfd.uni.geno
plot_ordination(F_HFD, F_weighted_uni_ord.HFD, type = "samples", color = "Genotype") + theme_bw(base_size = 14) +
theme(legend.position = "top") +
stat_ellipse() -> fig.hfd.w.uni.geno
fig.hfd.w.uni.geno
# Graphing out of curiosity
plot_ordination(F_HFD, F_uni_ord.HFD, type = "samples", color = "Surgery") + theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") + theme(legend.position = "top") + stat_ellipse()
plot_ordination(F_HFD, F_weighted_uni_ord.HFD, type = "samples", color = "Surgery") + theme_bw(base_size = 14) +
scale_color_brewer(palette="Dark2") + theme(legend.position = "top") + stat_ellipse()
plot_grid(
fig.hfd.uni.geno,
fig.hfd.w.uni.geno,
labels = "AUTO", ncol = 2, label_size = 20
) -> fig.hfd.summary
fig.hfd.summary
ragg::agg_jpeg(file.path(fig.dir, "fig_hfd_beta.jpeg"), width = 12, height = 5, units = "in", res = 600)
fig.hfd.summary
dev.off()
## png
## 2
sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] cowplot_1.1.1 ragg_1.2.6
## [3] doBy_4.6.19 DESeq2_1.40.2
## [5] SummarizedExperiment_1.30.2 Biobase_2.60.0
## [7] MatrixGenerics_1.12.3 matrixStats_1.1.0
## [9] GenomicRanges_1.52.1 GenomeInfoDb_1.36.4
## [11] IRanges_2.34.1 S4Vectors_0.38.2
## [13] BiocGenerics_0.46.0 phytools_1.9-16
## [15] maps_3.4.1 picante_1.8.2
## [17] nlme_3.1-163 ape_5.7-1
## [19] plotly_4.10.3 DT_0.29
## [21] ggrepel_0.9.4 microbiome_1.22.0
## [23] colorBlindness_0.1.9 RColorBrewer_1.1-3
## [25] pals_1.8 microViz_0.11.0
## [27] ggpubr_0.6.0 rstatix_0.7.2
## [29] vegan_2.6-4 lattice_0.22-5
## [31] permute_0.9-7 lubridate_1.9.3
## [33] forcats_1.0.0 stringr_1.5.0
## [35] dplyr_1.1.3 purrr_1.0.2
## [37] readr_2.1.4 tidyr_1.3.0
## [39] tibble_3.2.1 ggplot2_3.4.4
## [41] tidyverse_2.0.0 fs_1.6.3
## [43] phyloseq_1.44.0 magrittr_2.0.3
##
## loaded via a namespace (and not attached):
## [1] rstudioapi_0.15.0 jsonlite_1.8.7 farver_2.1.1
## [4] rmarkdown_2.25 zlibbioc_1.46.0 vctrs_0.6.4
## [7] multtest_2.56.0 RCurl_1.98-1.13 S4Arrays_1.2.1
## [10] htmltools_0.5.6.1 plotrix_3.8-2 broom_1.0.5
## [13] Rhdf5lib_1.22.1 rhdf5_2.44.0 gridGraphics_0.5-1
## [16] sass_0.4.7 bslib_0.5.1 htmlwidgets_1.6.2
## [19] plyr_1.8.9 cachem_1.0.8 igraph_1.5.1
## [22] lifecycle_1.0.3 iterators_1.0.14 pkgconfig_2.0.3
## [25] Matrix_1.6-1.1 R6_2.5.1 fastmap_1.1.1
## [28] GenomeInfoDbData_1.2.10 digest_0.6.33 numDeriv_2016.8-1.1
## [31] colorspace_2.1-0 crosstalk_1.2.0 textshaping_0.3.7
## [34] labeling_0.4.3 clusterGeneration_1.3.8 fansi_1.0.5
## [37] timechange_0.2.0 httr_1.4.7 abind_1.4-5
## [40] mgcv_1.9-0 compiler_4.3.1 microbenchmark_1.4.10
## [43] bit64_4.0.5 withr_2.5.1 doParallel_1.0.17
## [46] backports_1.4.1 BiocParallel_1.34.2 optimParallel_1.0-2
## [49] carData_3.0-5 ggsignif_0.6.4 MASS_7.3-60
## [52] DelayedArray_0.26.7 scatterplot3d_0.3-44 biomformat_1.28.0
## [55] tools_4.3.1 glue_1.6.2 quadprog_1.5-8
## [58] rhdf5filters_1.12.1 grid_4.3.1 Rtsne_0.16
## [61] cluster_2.1.4 reshape2_1.4.4 ade4_1.7-22
## [64] generics_0.1.3 gtable_0.3.4 tzdb_0.4.0
## [67] data.table_1.14.8 hms_1.1.3 Deriv_4.1.3
## [70] car_3.1-2 utf8_1.2.4 XVector_0.40.0
## [73] foreach_1.5.2 pillar_1.9.0 vroom_1.6.4
## [76] splines_4.3.1 bit_4.0.5 survival_3.5-7
## [79] tidyselect_1.2.0 locfit_1.5-9.8 Biostrings_2.68.1
## [82] knitr_1.44 xfun_0.40 expm_0.999-7
## [85] stringi_1.7.12 lazyeval_0.2.2 yaml_2.3.7
## [88] evaluate_0.22 codetools_0.2-19 cli_3.6.1
## [91] systemfonts_1.0.5 munsell_0.5.0 jquerylib_0.1.4
## [94] dichromat_2.0-0.1 Rcpp_1.0.11 mapproj_1.2.11
## [97] coda_0.19-4 ellipsis_0.3.2 bitops_1.0-7
## [100] phangorn_2.11.1 viridisLite_0.4.2 scales_1.2.1
## [103] crayon_1.5.2 combinat_0.0-8 rlang_1.1.1
## [106] fastmatch_1.1-4 mnormt_2.1.1